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ABSTRACT 



Aims. We extend earlier models of turbulent dynamos with an upper, nearly force-free exterior to spherical geometry, 
and study how flux emerges from lower layers to the upper ones without being driven by magnetic buoyancy. We also 
study how this affects the possibility of plasmoid ejection. 

Methods. A spherical wedge is used that includes northern and southern hemispheres up to mid-latitudes and a certain 
range in longitude of the Sun. In radius, we cover both the region that corresponds to the convection zone in the Sun 
and the immediate exterior up to twice the radius of the Sun. Turbulence is driven with a helical forcing function in 
the interior, where the sign changes at the equator between the two hemispheres. 

Results. An oscillatory large-scale dynamo with equatorward migration is found to operate in the turbulence zone. 
Plasmoid ejections occur in regular intervals, similar to what is seen in earlier Cartesian models. These plasmoid 
ejections are tentatively associated with coronal mass ejections (CMEs). The magnetic helicity is found to change sign 
outside the turbulence zone, which is in agreement with recent findings for the solar wind. 
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1. Introduction 

Observations show that the Sun shed s mass through 
twisted magnetic flux configurations (|Demoulin et all 



the turbulent convection. The dynamo cycle time can even 
be several hundred times the turnover time. In the corona, 
on the other hand, the typical time scale depends on the 
Alfven time, L/va, where L is the typical scale of mag- 
netic structures and va — -B/a/MoP is the Alfven speed for 
a given magnetic field strength B. Here, /io is the vacuum 
permeability. 

In a recent paper, IWarnecke fe Brandenburd (|2010T ) at- 
tempted a new approach of a unified treatment by combin- 
ing a dynamo-generated field in the convection zone with 
a nearly force-free coronal part, albeit in a local Cartesian 
geometry. In this paper, we go a step further by perform- 
ing direct numerical simulations (DNS) in spherical geom- 
etry. We also include density stratification due to grav- 
ity, but with a density contrast between the dynamo in- 
terior and the outer parts of the simulation domain that 
is much less (about 20) than in the Sun (around 14 orders 
of magnitude) . This low density contrast is achieved by us- 
ing an isothermal configuration with constant sound speed 
c s . Hence, the average density depends only on the gravita- 
tional potential and is given by lnp(r) w GM/rc 2 , where G 
is Newton's gravitational constant, M is the central mass, 
and r is the distance from the center. As convection is not 
possible in such an isothermal setup, we drive turbulence 
by an imposed helical forcing that vanishes outside the con- 
vection zone. This also helps achieving a strong large-scale 
magnetic field. The helicity of the forcing is negative (pos- 
itive) in the northern (southern) hemisphere and smoothly 
changes sign across the equator. Such a forcing gives rise to 
an a 2 dynamo with periodic o scillations and equa torward 

migration of magnetic activity (jMitra et all l2010ari . We ig- 

1 http://sohomre.nascom.nasa.gov/bestof soho/Movies/10thX«»^ia^^afl^gotation, SO there is no systematic shearing 
andhttp://www.youtube.com/watch?v=CvRj6Uykois&feature=pl^eiatatibelddeJi'ie only twisting comes then from the same 



2002) . Remarkable examples of such helical ejections can be 
seen in the movies produced by the SOHO and SDO mis- 
sionj]. Such events may be i mport ant for the solar dynamo 
(Blackm an fc Brandenburel 120031 ). They are generally re- 
ferred to as coronal mass ejections (CMEs). Conventionally, 
CMEs are modeled by adopting a given distribution of 
magnetic flux at the solar surface and letting it evolve by 
shearing an d twisting the magnetic field at its footpoints a t 
the surface (jAntiochos et all ll999t iTorok fc Klieml . l2003h . 
This approach is also used to model coronal he ating 
(|Gudiksen fe Nordlundl . [2001 iBingert fe Peterl . l20ll . The 
success of this method depends crucially on the ability to 
synthesize the velocity and magnetic field patterns at the 
surface. Of course, ultimately such velocity and magnetic 
field patterns must come from a realistic simulation of the 
Sun's convection zone, where the field is generated by dy- 
namo action. In other words, we need a unified treatment 
of the convection zone and the CMEs. The difficulty here 
is the large range of time scales, from the 11-year dynamo 
cycle to the time scales of hours and even minutes on which 
CMEs develop. Such a large range of time scales is related 
to the strong density stratification in the Sun, as can be 
seen from the following argument. In the bulk of the con- 
vection zone, the dynamo is controlled by rather slow mo- 
tions with turnover times of days and months. The typical 
velocity depends on the convective flux via F conv ps pu^msj 
where ~p is the mean density and « rms is the rms velocity of 
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motions that also sustain the dynamo-generated magnetic 
field. Note that in our model the mechanism of transport of 
the magnetic field to the surface is not magnetic buoyancy. 
Instead, we expect that, twisted magnetic fields will expel 
themselves to the outer regions by the Lorentz force. 

Our aim in this paper is not to provide a model as 
close to reality as possible, but to show that it is possi- 
ble to capture the phenomenon of CMEs (or, more gener- 
ally, plasmoid ejections) within a minimalistic model that 
treats the convection zone and the outer parts of the Sun 
in a self-consistent manner. That is, the magnetic field in 
the convection zone is dynamically generated by dynamo 
action and the motions are not prescribed by hand, but 
they emerge as a solution of the momentum equation and 
include magnetic backreaction from the Lorentz force. 

Given that gravity decreases with radius, there is in 
principle the possi bility of a radial win d with a critical point 
at r* = GM/2c 2 (IChoudhuril . Il998h . However, as we use 
stress-free boundary conditions with no mass flux in the 
radial direction, no such wind can be generated in our sim- 
ulations. Nevertheless, we observe radially outward prop- 
agation of helical magnetic field structures without mass 
flux. Furthermore, our results for the flux of magnetic he- 
licity compare we ll with recent measureme nts of the same 
in the solar wind (|Brandcu burg et all 1 2 llh . Our approach 
might therefore provide new insights not only for CMEs and 
dynamo theory, but also for solar wind turbulence. 



2. The model 

We use spherical polar coordina tes, (r,9,cj)). As in earlier 
work of lMitra et all (|2009L 120104 , our simulation domain is 
a spherical wedge. The results of lMitra et al.l (120091) for such 
a wedge are consistent with those o f lLivermore et al. (2010) 
for a full spherical shell, both of which ignored the effects of 
an equa tor, which was included in the work of lMitra et al.l 
( 2010a). Our model is a bi-layer in the radial direction. The 
inner layer is forced with random helical forcing functions 
which have different signs of helicity in the two hemispheres. 
This models the helical aspects of convection in the Sun. 
We shall often call the inner layer "turbulence zone" . The 
radius separating the two layers corresponds to the solar 
radius, r — R. This length scale is used as our unit of 
length. The inner layer models some aspects of the convec- 
tion zone (0.7R < r < R) without however having any real 
convection, and the outer layer (R < r < 2R) models the 
solar corona. We consider the range tt/3 < 9 < 27r/3, corre- 
sponding to ±30° latitude, and < (j) < 0.3, corresponding 
to a longitudinal extent of 17°. Here, is the polar angle 
and 4> the azimuth. At the solar surface at R — 700 Mm, 
this would correspond to an area of about 730 x 210 Mm 2 , 
which could encompass several active regions in the Sun. 
In our model the momentum equation is 



DU 

where F v 



-Vh 



J x B/p + F 



for 



(1) 



kinematic viscosity, Si 



p : V ■ (2pvS) is the viscous force, v is the 



j(Z7 i; j + Uj-i) - |%V ■ U is 
the traceless rate-of-strain tensor, semicolons denote co- 
variant differentiation, h — c 2 In p is the specific pseudo- 
enthalpy, c s = const is the isothermal sound speed, and 
g = —GMr/r 3 is the gravitational acceleration. We choose 
GM/Rc 2 = 3, so r* = 1.5R lies within our domain. This 



value is rather close to the surface and would lead to signif- 
icant mass loss if there was a wind, but this is suppressed 
by using impenetrative outer boundaries. 

The forcing function Ff or is given as the product of two 
parts, 



F for (r, 0, 0, i) =@ w (r- R) f(r, 9, <j>, t; 



cost 



(2) 



where Q w (r) = -^[1 — erf(r/w)] is a profile function con- 
necting the two layers and w is the width of the transi- 
tion at the border between the two layers (r = R). In 
other words, we choose the external force to be zero in 
the outer layer, r > R. The function / consists of ran- 
dom plane helical transverse waves with relative helicity 
a = (f ■ V X f)/k[f 2 and wavenumbers that lie in a band 
around an average forcing wavenumber of kfR ss 63. This 
value should be compared with the normalized wavenum- 
ber kiR, corresponding to the thickness of the shell AR, 
which yields k\R = 2irR/AR w 21, so the effective scale 
separation ratio, fcf/fci, is about 3. 

In Equation © the last argument of f(r 7 9,(f> 7 t;a) de- 
notes a parametric dependence on the helicity which is here 
chosen to be a = — cos 9 such that the kinetic helicity 
of the turbulence is negative in the northern hemisphere 
and positive in the southern. Specifically, / is given by 
(|Haueen et all I2003D 



f{x, t) = A f NBe{f k{t) exp[ifc(i) • x + i^(t)]}, 



(3) 



where Af is a nondimensional forcing amplitude, and x is 
the position vector. The wavevector k{t) and the random 
phase — 7r < cf)(t) < 7r change at every time step, so f(x, t) is 
(^-correlated in time. The normalization factor N is chosen 
on dimensional grounds to be N = c s {\k\c s / St) 1 / 2 . At each 
timestep we select randomly one of many possible wavevec- 
tors in a certain range around a given forcing wavenumber. 
The average wavenumber is referred to as k{. We ignore 
curvature effects in the expression for the forcing function 
and thus force the system with what would correspond to 
transverse helical waves in a Cartesian coordinate system, 
i.e., 



/ fe = R-/L nohcl) with R , ; *•< ; <7i -f \ 



(4) 



where — 1 < a < 1 is the helicity parameter of the forcing 
function, 



(nohcl) 

fc ~ 



(k x e) l^k 2 - (fc • e) 2 , 



(5) 



is a non-helical forcing function, and e is an arbitrary unit 
vector not aligned with fc; note that |/fc| 2 = 1. 

The pseudo-enthalpy term in Equation ([1]) emerges from 
the fact that for an isothermal equation of state the pressure 
is given by p = c 2 p, so the pressure gradient force is given 
by p~ 1 \ ! p = c 2 Vlnp = V/i. The continuity equation is 
then written in terms of h as 



Bh 



U. 



(G) 



Equations (dJ and © are solved together with the uncurled 
induction e quation for the ve c tor p otential A in the resis- 
tive gauge (|Candelaresi et all [201lh . 



dA 
~dt 



U xB + v W 2 A, 



(7) 
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where r\ is the (constant) magnetic diffusivity, so the 
magnetic field is given by B = V x A and thus 
obeys V ■ B — at all times. The gauge can in 
principle become important when calculating the mag- 
netic helicity density A ■ B, although the part result- 
ing from the small-s cale fields is expected to be inde- 
pendent of the gauge dSubramanian fe Brandenburg! [2006; 
iHubbard fc Brandenburg) . |201Q[) . while that of the large- 
scale fields is not. 

Our wedge is periodic in the azimuthal direction. For the 
velocity, we use stress-free boundary conditions on all other 
boundaries. For the magnetic field we employ vertical field 
conditions on r = 2R and perfect conductor conditions on 
both r — 0.7 R and the two 6 boundaries. Time is expressed 
in units of r = (u rms fcf) , which is the eddy turnover time 
in the turbulence zone, and it rms is the rms velocity in r < 
R. Density is given in units of the mean density in the 
turbulence zone, po — ~p. The magnetic field is expressed 
in units of the mean equipartition value, B eq , defined via 
B 2 q = popu 2 . We use a magnetic diffusivity that is constant 
in space and time and its value is given in terms of the 
magnetic Reynolds number, defined as 



Rm = u rms /r]k { 



(8) 



where k{ is the characteristic scale of the external force, 
defined above. This also turns out to be the energy con- 
taining scale of the fluid. In the following analysis, we 
use <j> averages, defined as F(r,6,t)=f F(r,6,(j>,t) d<f)/2w. 
Occasionally we also use time averages denoted by (.) t . We 
perform DNS with the Pencil Code0, which is a modu- 
lar high-order code (sixth order in space and third-order in 
time, by default) for solving a large range of partial differ- 
ential equations, including the ones relevant in the present 
context. 



3. Results 

3,1. Dynamo in the turbulence zone 

We start our DNS with seed magnetic field everywhere in 
the domain. Owing to the helical forcing in the turbulent 
layer, a large-scale magnetic field is produced by dynamo 
action. The dynamo is cyclic with equatorward migration 
of magnetic fields. This dynamo was studied by DNS in 
iMitra et al.l (|2010aF) and has been interpreted as an a 2 dy- 
namo. The possibility of oscillating a 2 dynamos was known 
sinc e the early papers of iBarvshnikova fc Shukurovl (|1987( ) 
and lRadler fc Brauerl (|l987h . who showed that a necessary 
condition for oscillations is that the a effect must change 
sign in the domain. 

The dynamo first grows exponentially and then satu- 
rates after around 300 turnover times, see Figure [T] After 
saturation the dynamo produces a large-scale magnetic 
field with opposite polarities in the northern and south- 
ern hemispheres. In Figure [5] we plot the radial magnetic 
field at the surface of the dynamo region at r = R. This 
layer would correspond to the solar photosphere if we had 
a more realistic solar model, which would include higher 
stratification as well as cooling and radiation effects. The 
six wedges represent different times and show clearly an 
equatorward migration of the radial magnetic field with 
polarity reversal every half cycle. The other components 
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t 1.5 



1000 2000 3000 4000 

t/T 



Fig. 1. Initial exponential growth and subsequent satura- 
tion behavior of the magnetic field in the interior for forced 
turbulence with dynamo action for Run A. The magnetic 
field strength is oscillating with twice the dynamo frequency 

*2uJ cvc . 



t/T=3028.00 t/r=3042.00 t/Y=3057.00 t/r=3072.00 t/T=3086.0G t/V=3l01.00 




http : //pencil-code . googlecode . com 



Fig. 2. Equatorward migration, as seen in visualizations of 
B r for Run D at r — R over a horizontal extent A8 = 58° 
and A0 = 17°. 

of the magnetic field (not plotted) also shows the same 
behavior. Comparing the first (t/r — 3028) and the last 
(t/r — 3101) panel, the polarity has changed sign in a time 
interval At/r ~ 100. The oscillatory and migratory prop- 
erties of the dynamo is also seen in the butterfly diagram 
of Figure [3] for (B$) r and (B r ) r . In Figure [T] one can also 
verify that the oscillation period is around 200 turnover 
times, corresponding to a non-dimensional dynamo cycle 
frequency of rw cyc = 0.032 and the field strength in the 
turbulent layer varies between 1.2 and 1.6 of the equipar- 
tition field strength. This value of th e cycle frequency is 
roughly consistent with an estimate of IMitra et al.l (|2010b|) 
that uj cyc = 0.5?7tfcinj where fc m is the relevant w avenumber 
of th e mean field. Using <q% ~ ?7to = «rms/3fcf (|Sur et all 
[2008f ). we find ™ cyc w 0.2(fc m /A:f) 2 w 0.02, where we have 
assumed fc m ~ 27r/0 3i? » 2 0fci and fcf 60/ci. The esti- 
mate of IMitra et al.l (|2010bl) applies to perfectly conduct- 
ing outer boundary conditions, which might explain the 
remaining discrepancy. 

A summary of all runs is given in Table [TJ where 
the amplitudes of the magnetic field show a weak non- 
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Fig. 3. Periodic variation of (B,f,) r and (B r ) r in the turbu- 
lence zone. Dark blue stands for negative and light yellow 
for positive values. The dotted horizontal lines show the 
location of the equator at — ir/2. The magnetic field is 
normalized by the equipartition value. Taken from Run A. 



monotonous dependence on the magnetic Reynolds number 
R m . For larger values of R m , the magnetic field strength 
decreases slightly with increasing i? m , but it is weaker 
than in some earlier a 2 dyna mos with open boundaries 
(jBrandenburg &: Doblerl . l200lh . This could be due to two 
reasons. Firstly, our simulations are far from the asymptotic 
limit of large magnetic Reynold s numb ers, in which the re- 
sults of iBrandenburg fc Doblerl (|200lh are applicable. The 
maximum value R m is in our simulations approximately 15 
times the critical i? m . The second reason could be that we 
have expulsion of magnetic helicity from our domain which 



Table 1. Summary of runs discussed in this paper. 



Run 


At 


Rm 


P m 


r>2 / n2 
"rms/ -°cq 


TCJ CyC 


At/T 




A 


0.20 


1.5 


1 


1.2-2.7 


0.032 


100 


0.482 


B 


0.20 


5 


1 


1.5-3.5 


0.029 


110 


0.409 


C 


0.20 


9 


1 


2.1-5.5 


0.022 


130 


0.455 


D 


0.20 


18 


1 


2.0-5.0 


0.019 


140 


0.409 


Dl 


0.10 


11 


1 


2.0-5.0 


0.018 


140 


0.455 


D2 


0.15 


15 


1 


2.0-5.0 


0.016 


130 


0.482 


D3 


0.25 


20 


1 


1.0-3.0 


0.023 


130 


0.293 


E 


0.20 


22 


1 


1.5-4.5 


0.017 


220 


0.205 


F 


0.20 


28 


1 


1.2-4.2 


0.015 


280 


0.273 


G 


0.20 


44 


1 


1.7-3.5 


0.015 


285 


0.409 
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Fig. 4. Time evolution of the radial magnetic field B r 
(solid line) and the azimuthal magnetic field B^ (dotted 
line) in the dynamo region averaged over the radius r and 
azimuth at 0=±7°. To improve the statistics, we calculate 
the components of the magnetic field as the antisymmetric 
part in latitude, i.e., B, = (Bf - Bf) /2 for i = r,<f>. 



Notes. Af is the forcing amplitude defined in Equation 
Rm is the magnetic Reynolds number defined in Equation ([8]) 
and P m — v jr\ is the magnetic Prandtl number. u cyc — 2n/T cyc 
stands for the frequency of the oscillating dynamo, where T cyc 
is the cycle period. At/r is the typical interval of plasmoid ejec- 
tions, whose typical speed is V e j. For the runs Dl to D3, the 3.4. Density variations 



was not present in (jBrandenburg fc Doblerl I2001D . We find 
the peak of the R m dependency at R m = 9, corresponding 
to Run C. The dynamo cycle frequency shows a weak de- 
crease (by a factor of 1.5) as the magnetic Reynolds number 
increases (by a factor of 20). 

3.2. Phase relation between radial and azimuthal fields 

Although our dynamo model does not include important 
features of the Sun such as differential rotation, some com- 
parison may still be appropriate. For the Sun, one measures 
the mean radial field by averaging the line-of-sight magnetic 
field from synoptic magnetograms. The azimuthal field is 
not directly observed, but its sign can normally be read 
off by looking at the magnetic field orientation of sunspot 
pairs. Existing data suggest that ra dial and azimutha l fields 
are approximately in out-of-phase ()Yoshimural I1976T ). This 
is reasonably well reproduced by afi dynamos models, 
where the radial field lags behind the azimuthal one by 
However, in the present work, radial and 
azimuthal fields are approximately in phase with a phase 
difference of 0.37T inside the dynamo region; see Figure |U 
Future studies will include the near-surface shear layer, 
which has been suspected to pla y an importan t role in pro- 
ducing equatorward migration (Brandenburg], l2005j ). This 
would also help reproducing the observed phase relation. 



3.3. Relation between kinetic and magnetic energies 

Next we investigate the relation between the rms values 
of the magnetic field and the velocity. Both quantities are 
oscillating in time with a typical period of 200 turnover 
times. In Figure [3] we compare the time evolution of the 
magnetic field strength and the rms velocity. The magnetic 
field is calculated in the dynamo region and normalized to 
the thermal equipartition field strength. The phase differ- 
ence between the two is slightly less than tt within the dy- 
namo region. This basically shows that the magnetic field 
quenches the turbulence. 



forcing amplitude A{ is changed, while r\ and v have the same 
value as for run D. The rms velocities in the turbulence zone 
change accordingly and affect therefore the Reynolds number 
and the turnover times r. 



The density is stratified in radius and varies by over an 
order of magnitude. For all the runs listed in Table [1] the 
density fluctuates about the hydrostatic equilibrium value, 
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Fig. 5. Phase relations between the magnetic field and the 
velocity in the dynamo region. The magnetic field is plot- 
ted as £?rms) normalized with the equipartition field of the 
(=[iopCg) as a solid and black line. The 



sound speed, 

rms velocity, normalized by the sound speed c s , is plotted 
as a dashed red line, and has been smoothed over 5 neigh- 
boring data points to make it more legible. Taken from 
Run A. 
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T 1 1 1 1 - 




_l i i i l_ 
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r/R 



Fig. 6. Radial dependence of density overplotted at differ- 
ent times. In the inset is the linear behavior of the loga- 
rithmic density logp/po to the inverse of the radius shown. 
Taken from Run A. 



p m po exp(GM/rc 2 ). The relative fluctuations are of com- 
parable strength at all radii; see Figure [SI 



3.5. Magnetic field outside the turbulence zone 

The magnetic field averaged over the entire domain is 
more than 5 times smaller than in the turbulence zone. 
In Figure [7] we show that the magnetic field is concentrated 
to the turbulence zone and B 2 drops approximately expo- 
nentially with a scale height of about 0.23i? in the outer 
parts for r > R. The toroidal component of the magnetic 
field is dominant in the turbulence layer, but does not play 
a significant role in the outer part. By contrast, the radial 
field is weak in the inner parts and dominates in the outer. 




0.8 1.0 1.2 1.4 1.6 1.8 2.0 
r/R 



Fig. 7. Radial dependence of the mean squared magnetic 
field, -Brms (solid line), compared with those of B r (dot- 
ted), Be (dashed), and B$ (dash-dotted). All quantities 
are averaged over 13 dynamo cycles. The inset shows the 
same quantities in a logarithmic representation. Taken from 
Run A. 



Magnetic structures emerge through the surface and 
create field line concentrations that reconnect, separate, 
and rise to the outer boundary of the simulation domain. 
This dynamical evolution is seen in a sequence of field line 
images in Figure HI where field lines of the mean field are 
shown as contours of rsmQA^ and colors represent B$. 
Prior to a plasmoid ejection we see a convergence of an- 
tiparallel radial field lines, which then reconnect such that 
the newly reconnected field lines move away from the re- 
connection site. The actual reconnection seems to happen 
much faster than the subsequent ejection. 

In the outer layers, the magnetic field emerges as large 
structures tha t corr e late w ith reconnection events of mag- 
netic fields. In iRustl ([19941 ) such phenomena have been de- 
scribed as magnetic clouds. We find recurrent ejections of 
magnetic field lines with concentrations and reconnection 
events, but the occurrence of structures such as magnetic 
clouds does not happen completely regularly, i.e., these 
structured events would be difficult to predict. 

In Figure [9] we show a close-up of the magnetic field. 
A configuration resembling a reconnection event is clearly 
seen. Here, the contours represent magnetic field lines 
with solid and dashed lines denoting counter-clockwise and 
clockwise orientations, respectively. The solid antiparallel 
field lines reconnect around r = 0.9R and separate to the 
left and to the right. On the right-hand side, a plasmoid has 
formed, which is eventually ejected. This plasmoid appears 
as a CME-like ejection in the first panel of Figure [TTJ These 
plasmoids, as seen more clearly in Figures 151 and 1TU1 appear 
as a concentration of field lines that propagate outwards. 
The fact that reconnection happens predominantly in the 
upper parts of the turbulence zone su ggests that turbulence 
is nee ded to enable fast reconnection (|Lazarian fc Vishniad 
llQQflh . 

Additionally, we find reconnection as a result of the in- 
teraction between ejections. As plotted in Figure [TU1 the 
field lines of two subsequent events have the opposite field 
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Fig. 8. Time series of formation of plasmoid ejections in spherical coordinates. Contours of rsin^A^, are shown together 
with a color-scale representation of B^; dark blue stands for negative and light yellow for positive values. The contours 
of r sin OA^ correspond to field lines of B in the rO plane. The dashed horizontal lines show the location of the surface 
at r = R. Taken from Run D. 
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Fig. 9. Time series of a reconnection event in an X-point as 
a close-up view. Contours of r sin 9Aj, are shown together 
with a color-scale representation of B^\ dark blue stands for 
negative and light yellow for positive values. The contours 
of r sin OA^ correspond to field lines of B in the rO plane, 
where solid lines represent counter-clockwise magnetic field 
lines and the dash ones clockwise. The dashed vertical lines 
show the location of the surface at r = R. Taken from Run 
D. 



line direction, which can then interact in the outer layers. 
Comparison with the first panel of Figure 1 1 X I shows that 
the current helicity has a correlation with the separatrices 
of the two polarities of the field lines. We also find that in 
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Fig. 10. Magnetic field configuration at the time of a 
ejection. Contours of rsinOA^, are shown together with a 
color-scale representation of B^\ dark blue stands for neg- 
ative and light yellow for positive values. The contours of 
r sin BA$ correspond to field lines of B in the rO plane, 
where the solid lines represent counter-clockwise magnetic 
field lines and the dash ones clockwise. The dashed vertical 
lines show the location of the surface at r — R. Taken from 
Run D. 



the interaction region the field lines have high density and 
are more strongly concentrated. 

3.6. Current helicity 

The current helicity (J • B) is often used as a useful proxy 
for the magnetic helicity (A ■ B) at small scales, because, 
unlike magnetic helicity, it is gauge- indep endent. Current 
helicity has also been observed in the Sun (jSeeha fcr. 1990) 
and it has been obtained from mean-field dynamo models 
(Dik pati &: Gilmanl . l200l[ ). In the present paper we are par- 
ticularly interested in the current helicity outside the Sun. 
We normalize it by the r-dependent time-averaged mean 
squared field to compensate for the radial decrease of J ■ B. 
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In Figure [T2J we have also averaged in latitude from 20° 
to 28°. In the turbulence zone the sign of J ■ B/(B 2 ) t is 
the same as that of kinetic helicity which, in turn, has the 
same sign as the helicity of the external forcing, i.e. of a; 
see Figures [TT] and IT2"1 

However, to our surprise, above the surface, and sepa- 
rately for each hemisphere, the signs of current helicity tend 
to be opposite to those in the turbulence zone; see Figure [TT1 
for the panels of t/r = 1669 and f/r = 1740. To demon- 
strate that plasmoid ejections are recurrent phenomena, we 
look at the evolution of J ■ B/(B 2 ) t as a function of t and 
r. This is done in Figure [12] for Run A. It turns out that 
the speed of plasmoid ejecta is about 0.45 of the rms ve- 
locity of the turbulence in the interior region for Reynolds 
numbers up to 15 and about 0.3 up to 30. To compare with 
the Sun, we estimate the rms velocity of the turbulence in 
terms of the convective energy flux via F « pu^ ms . The den- 
sity of the corona is p COI « 10 _12 kgm -3 , so our estimate 
would suggest V C) « 0.3(F/p cor ) 1/3 « 1200 km/s, which 
is somewhat above the observed speeds of 400-1000 km/s. 
The time interval between subsequent ejections is around 
100 t for Run A. As seen from Table [TJ the ejection interval 
is independent of the forcing amplitude At and increases 
weakly with magnetic Reynolds number, but it seems to 
be still comparable to half the dynamo cycle period, i.e., 
At « T cyc /2. This means that plasmoid ejections happen 
about twice each cycle. It is therefore not clear whether 
such a result can be extrapolated to the real Sun. 

In our simulations, we find the ejections to have the 
shape of the characteristic three-part structure that is ob- 
served in real CMEs. This is particularly clear in Figure [TTT 
where the ejections seem to contain three different parts. 
In the center we find a ball-shaped structure consisting of 
one polarity of current helicity, where at the front of the 
ejection a bow of opposite polarity had formed. In between 
these two structures the current helicity is close to zero and 
appears as a cavity. These thr ee pa r ts (pr ominence, cav- 
ity, and front) are described by lLowl ()1996l ) for CMEs and 
are generally referred to as 'three-part structure'. The basic 
shape of the ejection is independent of the used forcing am- 
plitudes and the kinetic and magnetic Reynolds numbers. It 
should, however, be noted that the three-part structure of 
the ejections becomes clearer at magnetic higher Reynolds 
numbers (e.g., for Runs D and G compared with Run A, 
for example). In the Sun, the plasma is confined to loops 
of magnetic field with flows along field lines due to the low 
plasma beta in the solar corona. This is also seen in our sim- 
ulations displayed in Figure [TT] where the ejections follow 
field lines and appear to create loop-like structures. An ani- 
mation of the detailed time evolution of the CME-like struc- 
tures emerging recurrently into the solar corona is available 
in the on-line edition (see Figure lll|P l However, since our 
choice of boundary conditions does not allow mass flux at 
the outer boundary, no plasma can actually leave the do- 
main. The re current nature of the plasm o id eje ctions found 
here and in IWarnecke fc Brandenburg! (l2010f) is not yet 
well u nderstood. In contrast to IWarnecke fc Brandenburg! 
(|2010| ). where there are no strong oscillations present, here 
the ejections seem to correlate with the dynamo cycle. In 
each hemisphere of the turbulence zone a magnetic field is 
created with different polarity. After they have migrated to 

3 The movie is also available at 

http : //www. youtube . com/watch?v=aR-PgxQyP24 
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Fig. 12. Dependence of the dimensionless ratio 
fioR J ■ B I (B 2 ) t on time t/r and radius r in terms 
of the solar radius. The top panel shows a narrow band in 
6 in the northern hemisphere and the bottom one a thin 
band in the southern hemisphere. In both plots we have 
also averaged in latitude from 20° to 28°. Dark blue stands 
for negative and light yellow for positive values. The dotted 
horizontal lines show the location of the surface at r = R. 



the equator, they merge and pro duce an ejection. However, 
compa ring with the results of IWarnecke fc Brandenburg] 
(|2010D . which are similar to those in the present paper, 
the oscillation cannot be the only explanation for the re- 
currence of the ejections. As we have seen in Figure [T2l 
these events export magnetic helicity out of the domain. 
For the dynamo, on the other hand, magnetic helicity losses 
play a role only in the nonlinear stage. It is therefore con- 
ceivable that the regular occurrence of plasmoid ejections 
is connected with nonlinear relaxation oscillations rather 
than with the dynamo cycle which is essentially a lin- 
ear phenomenon. This is also su ggested by the results of 
IWarnecke fc Brandenburg! (|2010T) . where recurrent ejections 
occur without oscillations of the large-scale field. 

From Figures !TT1 and !T2l we conclude that in each hemi- 
sphere the sign of current helicity outside the turbulence 
zone is mostly opposite to that inside the turbulence zone. 
A stronger trend is shown in the cumulative mean of cur- 
rent helicity over time. This is shown in Figure 1131 where 
we plot the time evolution of the <f> averaged current helic- 
ity at r — 1.5 R and 28° latitude, which is a safe distance 
away from the outer r and 9 boundaries so as not to perturb 
our results, which should thus give a reasonable represen- 
tation of the outer layers. For the northern hemisphere the 
current helicity (solid black line) and the cumulative mean 
(solid red line) show positive values and for the southern 
hemispher e (dotted lines) neg a tive v alues. This agrees with 
results of iBrandenburg et al.l (]2009() . where the magnetic 
helicity of the field in the exterior has the opposite sign 
than in the interior. 

To investigate whether the sign of the current helicity is 
different in the turbulence zone and in the outer parts, we 
show in Figure [14] the azimuthally and time-averaged cur- 
rent helicity as a function of radius and colatitude. It turns 
out that, even though we have averaged the result over sev- 
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Fig. 11. Time series of coronal ejections in spherical coordinates. The normalized current helicity, {IqRJ ■ B/(B 2 ) t , is 
shown in a color-scale representation for different times; dark blue stands for negative and light yellow for positive values. 
The dashed horizontal lines show the location of the surface at r = R. Taken from Run D. The temporal evolution is 
shown in a movie available as online material. 




Fig. 13. Dependence of the dimensionless ratio 
HoR J ■ B/(B 2 ) t on time t/r at radius r — 1.5 i? and 
28° latitude. The solid line stands for the northern hemi- 
sphere and the dotted for the southern hemisphere. The red 
lines represent the cumulative mean for each hemisphere. 



eral thousand turnover times, the hemispheric sign rule of 
current helicity is still only approximately obeyed in the 
outer layers — even though it is nearly perfectly obeyed in 
the turbulence zone. Nevertheless, there remains substan- 
tial uncertainty, especially near the equator. This suggests 
that meaningful statements about magnetic and current he- 
licities in the solar wind can only be made after averaging 
over sufficiently long stretches of time. 

3. 7. Magnetic helicity fluxes 

In view of astrophysical dynamo theory it is important to 
understand the amount of magnetic helicity that can be 
exported from the system. Of particular interest here is the 
magnetic helicity associated with the small-scale magnetic 
field. Under the assumption of scale separation, this quan- 
tity is gauge-independent (|Subramanian fc Brandenburgj 
2006), so we can express it in any gauge. T his has been 
verifie d in simulations b oth with an equator (|Mitra et all 
I2010bf) and without (jHubbard fe Brandenbursd |2010D . 
Here, we compute the magnetic helicity flux associated with 



Fig. 14. Current helicity averaged over 3900 turnover 
times. Legend is the same as in Figure [TTJ Dark blue cor- 
responds to negative values, while the light yellow corre- 
sponds to positive value. Taken from Run D. 

the small-scale field by subtracting that of the azimuthally 
averaged field from that of the total field, i.e., 



exa = ExA — E x A, 



(9) 



where E — /lorjJ—UxB is the electric field. This is also the 
way how the magn etic helicity flux from the small-sc ale field 
was computed in iHubbard fc Brandenburg (|2010t ). where 
the magnetic helicity flux from the total and large-scale 
fields was found to be gauge-dependent, but that from the 
small-scale field was not. In Figure HH] we compare the flux 
of magnetic helicity of the small-scale field across the outer 
surfaces in the northern and southern hemispheres with 
that through the equator. It turns out that a major part of 
the flux goes through the equator. The part of the magnetic 
helicity flux that goes through the surface is about 20% of 
r) t B 2 q . However, the magnetic helicity flux should primarily 

depend on B rather than B eq . In the present simulations, 
where the dynamo works with a fully helical field, the two 
are comparable. This is not the case in the Sun, where 
the estimated f ield s trength is expected to be about 300 G 
(|Brandenburei I2009D . Thus, to compare with the Sun, a 
more reasonable guess for the magnetic helicity flux would 
be about 20% of rjtB 2 . Integrating this over one hemisphere 
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Fig. 15. Time evolution of the magnetic helicity flux of 
the large-scale field, smooth over two data points. Here, 
the mean of magnetic helicity flux out through the surface 
of the northern hemisphere (black) is shown, together with 
that through the southern hemisphere (dotted red), and the 
equator (dashed blue). 



and multiplying this with the 11 year cycle time, we find 
the total magnetic helicity loss to be 2irR 2 ri t B 2 T cycl which 
corresponds to 5 x 10 47 Mx 2 if we use i) t = 3x 10 12 cm 2 s _1 . 
This value exceeds the estimated upper limi t for the solar 
dynam o of about 10 46 Mx 2 per cycl e given bv | Brandenburgj 
(|2009fl . However the estimate by iBrandenburd (|2009t ) is 
based on an afl dynamo model with a effect and shear 
that yield a period comparable with the 11 year period of 
the Sun. Therefore, the discrepancy with the present model, 
where shear plays no role, should not be surprising. Instead, 
it tells us that a dynamo without shear has to shed even 
more magnetic helicity than one with shear. 

The magnetic helicity flux of the large-scale field may 
not a priori be gauge-invariant. However, the system is in 
a statistically steady state and, in addition, the magnetic 
helicity integrated over each cycle is constant. In that case 
the divergence of the magnetic helicity flux is also gauge- 
invariant. Furthermore, the shell-integrated magnetic he- 
licity cannot have a rotational component and is therefore 
uniquely defined. In Figure[15]we plot this flux and see that 
its maxima tend to occur about 50 turnover times after 
magnetic field maxima; see Figures H] and [5j The helicity 
flux of the small-scale field does not show a clear behav- 
ior. Since the ejections appear to be related to the mag- 
netic field strength in this way, one might conclude that 
the magnetic helicity flux of the large-scale field is trans- 
ported through these ejections. This result is somewhat un- 
expected and deserves to be reexamined more thoroughly 
in future simulations where cycle and ejection frequencies 
are clearly different from each other. 

Next, let us look at the magnetic helicity flux of the 
small-scale field. On earlier occasions, iMitra et all (l2010bh 
and lHubbard k, Brandenburg! (|2010D have been able to de- 
scribe the resulting magnetic_helicity flux_by a Fickian dif- 
fusion ansatz of the form Ff = — KhV/if, where Kh/vto 
was found to be 0.3 and 0.1, respectively. In Figure [T7] 
we show that the present data allow a similar representa- 
tion, although the uncertainty is large. It turns out that 
Kh/vto is about 3, suggesting thus that turbulent mag- 



Fig. 16. Cumulative mean of the time evolution of the 
magnetic helicity flux of the small-scale field, F[ = e x a, 
normalized by ?7t-B 2 q , where rj t ~ t^o = «rms/3fcf was de- 
fined in Section 13.11 Here, the mean of magnetic helicity 
flux out through the surface of the northern hemisphere 
(black) is shown, together with that through the southern 
hemisphere (dotted red), and the equator (dashed blue). 
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Fig. 17. Dependence of the latitudinal component of the 
— f 

magnetic helicity flux, F s , compared with the latitudinal 
gradient of the magnetic helicity density of the small-scale 
field, Ve»/if, at r/R — 0.85. The latter agrees with the for- 
mer if it is multiplied by an effective diffusion coefficient for 
magnetic helicity of Kt ~ 3?7to- 



netic helicity exchange across the equator can be rather 
efficient. Such an efficient transport of magnetic helicity 
out of the dynamo region is known to be beneficial for 
the dynamo in that it alleviates catastrophic quenching 
(|Blackman fc Brandenburel 120031 ). In this sense, the inclu- 
sion of CME-like phenomena is not only interesting in its 
own right, but it has important beneficial consequences for 
the dynamo itself in that it models a more realistic outer 
boundary condition. 

3.8. Comparison with solar wind data 

Our results suggest a reversal of the sign of magnetic he- 
licity between the inner and outer parts of the computa- 
tional domain. This is in fact in agreement with recent 
attempts to measure magnetic helicity in the solar wind 
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Fig. 18. Hclicity in the northern outer atmosphere. The 
values are written out at the point, r = 1.5 R, 90° — 9 — 17°, 
and = 9°. Top panel: Phase relation between the toroidal 
and poloidal Bg field, plotted over time t/r. Bottom 
panel: Helicity H(k) is plotted over normalized wave num- 
ber kR. The helicity is calculated with the Taylor hypoth- 
esis using the Fourier transformation of the poloidal and 
toroidal field. 



(|Brandenburg et all . 1201 lh . They used the Taylor hypoth- 
esis to relate temporal fluctuations of the magnetic field to 
spatial variations by using the fact that the turbulence is 
swept past the space craft with the mean solar wind. This 
idea can in principle also be applied to the present simula- 
tions, provided we use the obtained mean ejection speed V^j 
(see Table [T]) for translating temporal variations (in t) into 
spatial ones (in r) via r = rq— V e jt. Under the assumption of 
homogeneity, one can then estimate the magnetic helicity 
spectrum as H(k) = Alm(BgB^)/k; see iMatthaeus et al.l 
<|!982h and Eq. (9) of lBrandenburg et al.l (|201lh . Here, hats 
indicate Fourier transforms and the asterisk denotes com- 
plex conjugation. 

In Figures [TB] and \W\ we show the results for the north- 
ern and southern hemispheres, as well as time series of the 
two relevant components Bg and B^. The resulting mag- 
netic helicity spectra, normalized by 2yLQEyi/k, where Em 
is the magnetic energy spectrum, give a quantity that is 
between —1 and +1. Note that the time traces are gov- 
erned by a low frequency component of fairly large ampli- 
tude. In addition, there are also other components of higher 
frequency, but they are harder to see. The results suggest 
positive magnetic helicity in the north and negative in the 
south, which would be indicative of the helicities of the 
solar wind at smaller length scale. It also agrees with the 
current helicities determined using explicit evalu ation in 
real s pace. On the other hand, the Parker spiral (jParkerL 
Il958f ) might be responsible f or the magnetic helicity at large 
scales (|Bieber et all Il987allbh . 



South 

0.010 c 




-o.oio r 

260 280 300 320 340 360 380 

t/r 




Fig. 19. Helicity in the southern outer atmosphere. The 
values are written out at the point, r = 1.5 R, 90° — 9 = 
— 17° and <j> — 8.6°. Top panel: Phase relation between 
the toroidal B^ and poloidal Bg field, plotted over time 
t/r. Bottom panel: Helicity H(k) is plotted over normal- 
ized wave number kR. The helicity is calculated with the 
Taylor hypothesis using the Fourier transformation of the 
poloidal and toroidal field. 

4. Conclusions 

In the present work we have demonstrated that CME- 
like phenomena are ubiquitous in simulations that include 
both a helicity-driven dynamo and a nearly force-free exte- 
rior above it. This was firs t show n in Cartesian geometry 
([Warnecke fc Brandenburel I2010T ) and is now also verified 
for spherical geometry. A feature common to both mod- 
els is that the helical driving is confined to what we call 
the turbulence zone, which would correspond to the con- 
vection zone in the Sun. In contrast to the earlier work, 
we have now used a helical forcing for which the kinetic 
helicity changes sign across the equator. This makes the 
dynamo oscill atory and displays e quatorward migration of 
magnetic field iMitra et al.l (|2010al) . More importantly, un- 
like our earlier work where the gas pressure was neglected 
in the outer parts, it is fully retained here, because it does 
automatically become small away from the surface due to 
the effect of gravity that is here included too, but was ne- 
glected in the ear lier Cartesian model. The solution s shown 
here and those of lWarnecke fc Brandenb urg (2010) demon- 
strate that this new approach of combining a self-consistent 
dynamo with a corona-like exterior is a viable one and can 
model successfully features that are similar to those in the 
Sun. However, our model is still not sophisticated enough 
for direct quantitative comparisons. 

Of particular interest is the sign change of magnetic 
and current helicities with radius. Although similar be- 
havior has also be e n see n in other Cartesian models of 
iBrandenburg et ahl ((20091. its relevance for the Sun was 
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unknown until evidence for similar sign properties emerged 
from solar wind data (jBrandenburg et all 1201 if) . In the 
present case we were also able to corroborate similar find- 
ings by using the Taylor hypothesis based on the plasmoid 
ejection speed. It is remarkable that this appears to be suf- 
ficient for relating spatial and temporal fluctuations to each 
other. 

There are many ways in which the present model can be 
extended and made more realistic. On the one hand, the as- 
sumption of isothermal stratification could be relaxed and 
the increase of temperature in the corona together with 
the solar wind could be modeled in a reasonably realis- 
tic way. On the other hand, the dynamo model could be 
modified to include the effects of convection and of latitu- 
dinal differential rotation. Among other thi ngs, diff e rentia l 
rotation would lead to the Parker spiral (jParkerL [l958), 
which is known to pr oduce magnetic helicity of its own 
(jBieber et all ll987 a.b N l . It would then be interesting to see 
how this affects the magnetic helicity distribution seen in 
the present model. 
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